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ABSTRACT 



Aims. We investigate the possible nonlinear variability properties of the black hole X-ray nova 4U 1543-47 with a dual 
goal: 1) to complement the temporal studies based on linear techniques, and 2) to search for signs of (deterministic 
and stochastic) nonlinearity in Galactic black hole (GBH) light curves. The proposed analysis may provide additional 
model-independent constraints to shed light on black hole systems and may strengthen the unification between GBHs 
and Active Galactic Nuclei (AGN). 

Methods. First, we apply the weighted scaling index method ( WSIM) to characterize the X-ray variability properties of 
4U 1543-47 in different spectral states during the 2002 outburst. Second, we use surrogate data to investigate whether 
the variability is nonlinear in any of the different spectral states. 

Results. The main findings from our nonlinear analysis can be summarized as follows: 1) The mean weighted scaling 
index (a) appears to be able to parametrize uniquely the temporal variability properties of GBHs. The 3 different 
spectral states of the 2002 outburst of 4U 1543-47 are characterized by different and well constrained values of (a) 
satisfying the following relationship: (q)vhs < (q)ls < (q)hs- 2) The search for nonlinearity reveals that the 
variability is linear in all light curves with the notable exception of the very high state (VHS). 

Conclusions. Our results imply that we can use the WSIM to assign a single number, namely the mean weighted scaling 
index (a), to a light curve, and in this way discriminate among the different spectral states of a source. The detection 
of nonlinearity in the VHS, that is characterized by the presence of most prominent QPOs, suggests that intrinsically 
linear models which have been proposed to account for the low frequency QPOs in GBHs are probably ruled out. Finally, 
the fact that WSIM results are scarcely affected by the noise level and length of the light curve, naturally suggests 
an application to AGN variability with the possibility of a direct comparison with GBHs. However, before deriving 
more general conclusions, it is first necessary to carry out a systematic nonlinear analysis on several GBHs in different 
spectral states in order to assess whether the results obtained for 4U 1543-47 can be considered as representative for 
the entire class of GBHs. 

Key words. Methods: data analysis - X-rays: binaries - X-rays: individuals: 4U 1543-47 



1. Introduction 

In recent years, several studies have demonstrated the im- 
portance of X-ray temporal and spectral studies of black 
hole systems. Because of their closeness and brightness, the 
physical conditions of Galactic black holes (GBHs) are bet- 
ter known than those of supermassive black holes in active 
galactic nuclei (AGN), and in principle can be used to in- 
fer information about their scaled-up extragalactic analogs. 
For example, is it now well accepted that GBHs undergo a 
continuous spectral evolution, switching between two main 
states: the "low/hard" (hereafter LS) and the "high/soft" 
(HS) passing through less well-defined and short-lived "in- 
termediate states" , sometimes called steep power law, SPL, 
state or very high state, VHS, if occurring at high luminos- 
ity (see McClintock & Remillard 2006, and Done et al. 2007 
for recent comprehensive reviews on GBHs). 

Despite a substantial advance in this field, however, sev- 
eral questions yet remain unanswered. For example, it is 



currently strongly debated whether the LS is character- 
ized by a truncated disk or not (e.g., Miller et al. 2006a, b; 
Gierlinski et al. 2008; D'Angelo et al. 2008). It is still un- 
known if the jet plays an important role in the X-ray range 
during the LS (e.g., Markoff et al. 2001; Zdziarski et al. 
2003), and whether there is a physical difference between 
the LS and the quiescent state (e.g., Tomsick et al. 2004; 
Corbel et al. 2006). It is also still unclear what is the ori- 
gin of QPOs, which only appear in specific spectral states. 
Finally, we still have a poor knowledge of the physical con- 
ditions of the accretion flow in the VHS, which is often as- 
sociated with the most powerful relativistic ejections (e.g., 
Fender et al. 2004). 

Several different models, mostly driven by X-ray spec- 
tral results, have been proposed to explain the aforemen- 
tioned open questions. However, due to the transient nature 
and short duration of these phenomena, the spectral infor- 
mation alone is unable to discriminate between competing 
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In order to break this degeneracy, it is therefore of cru- 
cial importance to complement the spectral information 
with additional constraints from the temporal analysis. In 
this framework, the use of the power spectral density (PSD) 
functions has proved to be very successful: the combined 
temporal and spectral information has led to a generally 
accepted scenario where the spectral evolution of GBHs is 
mostly driven by variations of the accretion rate m, which 
in turn lead to changes in the interplay between accre- 
tion disk, Comptonizing corona, and a relativistic jet (e.g., 
McClintock & Remillard 2006). 

The success of the PSD in GBH and AGN studies (see, 
e.g., Klein- Wolt & van der Klis 2008; McHardy et al. 2006), 
emphasizes the importance temporal studies and the crucial 
role that they may play in breaking the spectral degeneracy. 
It is however important to explore also alternative temporal 
methods since the PSD, as any other timing techniques, 
cannot exhaustively characterize any non-trivial variable 
system. For example, since the PSD is sensitive only to the 
first 2 moments of the probability distribution, it cannot 
provide a complete description of non-Gaussian processes. 
Similarly, since the PSD is an intrinsically linear technique 
it is not adequate to characterize systems with non-linear 
variability. 

The main goal of this work is to investigate alternative 
temporal methods, that are complementary to the PSD. 
More specifically, in the first part of this paper (§3, and §4), 
we carry out a nonlinear analysis of the variability proper- 
ties of the black hole X-ray nova 4U 1543-47 (whose data 
description is provided in §2) by using the weighted scal- 
ing index method (WSIM). The scaling index method (e.g., 
Atmanspacher et al. 1989) has been employed in a number 
of different fields because of its ability to discern underlying 
structure in noisy data. For example, it has been success- 
fully used in medical science (see Morfill & Bunk 2001 for a 
brief review), in image analysis (e.g., Rath & Morfill 1997; 
Brinkmann et al. 1999), in plasma physics (e.g., Ivlev et al. 
2008; Sutterlin et al. 2009), and in cosmology (Rath et al. 
2002, 2007, 2009). Recently, we have applied this method 
to well-sampled AGN light curves to look for signs of non- 
stationarity (Gliozzi et al. 2002; 2006). 

In this work, our primary goal is to investigate whether 
we can parametrize the variability properties of a GBH with 
a single number, similarly to what is currently done in the 
spectral analysis where the energy spectra of accreting ob- 
jects are usually characterized by one number, namely the 
slope of the power-law model. To this aim, we apply the 
WSIM to 4U 1543-47, a BH system for which the corre- 
spondence between spectral states and individual observa- 
tions is well defined during its outburst in June 2002. If 
the results from this analysis are encouraging, we plan to 
carry out a similar systematic analysis on a sample of BH 
systems, to investigate whether a common pattern emerges. 

The second part of this works (§5) deals with the search 
for nonlinearity in the light curves of 4U 1543-47, in all its 
spectral states. We use a new method to produce reliable 
surrogate data and the nonlinear prediction error (NLPE; 
Sugihara & May 1990) test to search for any signs of non- 
linearity. As discussed below, this statistic is able to detect 
any non-linear behavior in the light curves, irrespective of 
whether the non-linearity is "deterministic" or "stochastic" 
in nature. 

It is worth noting that, unlike previous nonlinear studies 



search for low-dimensional chaotic signatures in GBHs, but 
to characterize the global variability properties of BHs in a 
simple and scalable way. We stress that this kind of analysis 
is not in contrast with the "standard" linear analysis, whose 
contribution is of crucial importance in guiding our work, 
but rather it complements it by exploring aspects of the 
variability that are not accessible to linear techniques. 

2. Data Description 

4U 1543-47 is a recurrent X-ray nova, with outbursts oc- 
curring every 10-12 years. Dynamical optical studies dur- 
ing quiescence yield a primary mass of Mi = 9.4 ± 2.0M Q 
strongly arguing for the presence of a BH (Park et al. 2004) . 

For our analysis we will use RXTE PCA data of 
4U 1543-47 during the outburst occurred between June 
17 and July 22, 2002, which corresponds to the interval 
52,442-52,477 in Modified Julian Date (MJD= Julian Date- 
2,400,000.5). The RXTE PCA covered the 2002 flare of 
4U 1543-47 with at least one observation per day with ex- 
posures ranging between ^800 s and ^3900 s. More than 
90% of the observations caught the source in high state 
(HS), whereas only a couple of observations cover the short- 
lived very high state (VHS) and the beginning of the low 
state (LS). 

A detailed analysis of the energy and power spectral 
densities was performed by Park et al. (2004). For our pur- 
pose, the relevant results of their work can be summarized 
as follows: 1) The source was caught by the PCA close to 
the outburst peak, when 4U 1543-47 was in a thermally 
dominated state HS. 2) For nearly the entire duration of 
the outburst, the source remained in the HS, which is tem- 
porally characterized by a featureless PSD (see Figs. 8a, c of 
Park et al. 2004). The only notable exceptions are: a rapid 
transition to the VHS around 52,459-52,460 MJD, whose 
PSD shows a very prominent QPO around 5-10 Hz (Fig. 
8b of Park et al. 2004) , and a transition to the LS toward 
the end of the observation showing the typical band-limited 
noise PSD (Fig. 8d of Park et al. 2004). 3) The energy spec- 
tral analysis indicates that an acceptable parametrization 
of all spectra always requires a disk, a power-law compo- 
nent, and a Fe Ka line, whose relative contributions vary 
with time. Specifically, the disk flux closely follows the total 
count rate evolution during the outburst, whereas the flux 
associated with the power-law component shows a broad 
and prominent peak during the VHS, preceded by 2 less 
prominent peaks. 

All light curves were extracted in the 2-20 keV energy 
range following the standard RXTE procedure and were 
binned at a resolution of 0.100097656 s (for simplicity here- 
after we will use 0.1s when referring to the time bin), which 
is 205 times the original integration time. A more detailed 
description of the data reduction is provided by Reig et al. 
(2006). 

3. Weighted Scaling Index Method 

Since a detailed description of the scaling index method 
(SIM) has already been provided in our previous work, 
here we limit ourselves to summarizing the main steps in a 
simple and qualitative way and pointing out the main dif- 
ferences introduced by the weighted scaling index method 
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Fig. 1. Representative light curves during HS (top panel), 
VHS (middle panel), and LS (bottom panel). Time bin is 
0.1 s 



1. As all timing techniques, the very first step starts with a 
time series. Figure 1 shows three samples of light curves 
characterizing HS (top panel), VHS (middle panel), and 
LS (bottom panel) of 4U 1543-47 during the 2002 out- 
burst. For clarity reasons, 500 s time intervals (i.e., in- 
tervals with 5000 data points) are plotted, although 
the results described in §4.1 are obtained using 1000 
s (10,000 points) intervals; the time bin used for all 
light curves is 0.1 s. Before applying the SIM, the light 
curves are normalized, in the sense that the mean count 
rate is subtracted from each data point and the result- 
ing quantity is divided by the total standard deviation. 
In any variability analysis, particular attention should 
be paid to low count rate intervals due to the Poisson 
noise that may swamp the signal and hamper the analy- 
sis. Despite the large difference in count rate and hence 
in Poisson noise (during the 2002 outburst the average 
RXTE PCA count rate during individual observations 
ranges between ~18, 000-400 count/s), the SIM analy- 
sis appears to be largely unaffected by the actual count 
rate (see §4.2 for more details). 

2. As most of the methods of nonlinear dynamics applied 
to timing analysis, the SIM relies upon the phase space 
reconstruction, which is obtained via the time-delay 




Fig. 2. 3-dimensional phase space portraits for HS (top 
panel), VHS (middle panel), and LS (bottom panel), r =0.1 
s was used for the time-delay reconstruction. 



2006). In simple words, in the case of a 3-dimensional 
(3-D) phase space reconstruction, one constructs a set of 
3-D vectors by selecting triplets of data points from the 
original time series, which are separated in time by the 
time-delay r. More specifically, the second data-point 
(which represents the y-component in a 3-D vector) is 
separated from the first one (the x-component) by a 
time delay of r, whereas the third data-point (the z- 
component) is separated from the first one by 2 r, and 
from the second one by r. Figure 2 illustrates the phase 
space portraits in 3 different spectral states, i.e. the re- 
sults of a 3-D phase space reconstruction for the light 
curves shown in Figure 1 using a time delay of r = 0.1 s. 
Not surprisingly, just like the light curve plotted in the 
middle panel of Fig. 1 looks different from the other two, 
the phase space portrait of the VHS looks considerably 
different from those describing HS and LS. This reflects 
the fact that, unlike HS and LS, the VHS light curve is 
characterized by the presence of frequent peaks and dips 
in a quasi-periodic fashion, which increase the degree of 
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3. After the phase space reconstruction, the temporal 
properties of the original time series translate into topo- 
logical properties of the phase space portraits. In or- 
der to quantify these topological properties, a useful 
method is based on the measure of the crowding around 
each individual vector. This measurement is formally 
performed by computing the cumulative number func- 
tion, d{R) = n{j\dij < R} 7 which measures the num- 
ber of vectors j, whose distance dij from a vector i 
is smaller than R. Generally, when plotted versus R 
in log- log space, the function Ci{R) can be approxi- 
mated by a power law, Cj (R) oc R ai , for a wide range 
of values of R. The exponent, on = [log((7j(i?2)) — 
log(Ci(i?i))]/[log(i? 2 ) - log(-Ri)], which is the logarith- 
mic derivative of the cumulative number function, is the 
scaling index. In summary, for a light curve with N data 
points and an embedding space of dimension D this pro- 
cess yields N — D values of on. The temporal properties 
of the original light curve can then be studied either 
using the distribution of on or the mean value of this 
distribution (a), as explained in §4. 

4. The main difference between "normal" and weighted 
scaling index methods is that in the latter the cumu- 
lative number function is substituted by the weighted 
cumulative point distribution, where the weighting func- 
tion can be any diffcrentiable function (in our case a 
Gaussian function; see Rath et al. 2007 for a detailed ex- 
planation of the WSIM). The chief advantage of WSIM 
is twofold: 1) the logarithmic derivative (i.e., the scaling 
index) can be computed analytically instead of numer- 
ically; 2) the number of free parameters is reduced by 
1: since the logarithmic derivative is computed analyti- 
cally, we only need to define one value R at which it is 
computed, instead of R\ and i?2- 

To summarize, i) we start from normalized light curves, 
ii) transform them into phase space portraits via time-delay 
reconstruction, and iii) quantify their differences by com- 
puting the weighted scaling index. 

Importantly, the mean value (a) provides direct infor- 
mation on the nature of the variability process: for a purely 
random process (a) tends to the value of the dimension of 
the embedding space (i.e., the space used in the phase space 
reconstruction) , whereas for correlated (and deterministic) 
processes (a) is always smaller than the dimension, D, of 
the embedding space and, in the ideal case where the ran- 
dom component is completely negligible, the mean scaling 
index is independent of the embedding dimension. In other 
words, low values of (a) characterize correlated variability 
processes, whereas higher values correspond to variability 
properties with a higher degree of "randomness" . 

Before discussing the results of the WSIM, it is impor- 
tant to understand the role played by the 3 free parameters 
involved in this process and their impact on the results. The 
process leading to the phase space reconstruction requires 
to choose 2 parameters, the time delay r and the embed- 
ding dimension D, and the computation of the weighted 
scaling index requires an additional parameter, the radius 
R at which the logarithmic derivative is computed. 

A common choice for r is the characteristic timescale of 
the system, which can be determined with different meth- 
ods (e.g., autocorrelation function, PSD, or the so-called 
mutual information; Fraser & Swinney 1986). In our anal- 



4U 1543-47 show prominent QPOs in the VHS and LS, 
whereas as expected the PSD of the HS is featureless (see 
Park et al. 2004). On the other hand, there are no system- 
atic prescriptions for the choice of D and R. The latter 
likely depends on the typical distances (in the following we 
will use the Euclidean norm as measure of the distance) 
between vectors, which, in turn, depend on the choice of 
the embedding space dimension. 

In principle, the discriminating power of the SIM is 
enhanced by using high embedding dimensions. However, 
for our purposes, a relatively low embedding dimension is 
preferable, since we work with a limited number of points 
and since one of our goals is to compare the GBHs vari- 
ability results with those of AGN, which have light curves 
characterized by fewer data points than GBHs. Specifically, 
for our analysis we utilize D = 3. Once t and D are fixed, 
R is obtained in the following way: first, the temporal order 
of the data points in the original time series is randomized 
creating 10 sets of randomized data; second, the mean WSI 
is computed for randomized and original data using differ- 
ent values of R (specifically between 0.5 and 2.5); finally, 
the chosen radius (in our case R = 1.6) is the value that 
yields the larger difference between randomized data and 
original data, indicating that it is the value most sensitive 
to the temporal structure of the original data. 

Our nonlinear analysis of the variability properties of 
4U 1543-47 is therefore carried out using r = 0.1 s, D = 3, 
and R — 1.6. However, for the sake of completeness, we 
have performed an investigation of a broad parameter space 
encompassing r = 0.1 - 10 s, D = 2 - 4, R = 0.6 - 2.4. 
The results of this analysis, which demonstrates that our 
main findings are largely insensitive to the choice of these 
3 parameters, are reported in the Appendix. 



4. WSIM Results 

As explained before, the temporal properties of a variable 
system can be studied either via the distribution of the 
on values or simply through the mean value (a). In the 
following we elucidate these 2 approaches, by applying the 
WSIM first to the 3 representative light curves shown in 
Figure 1, and then to all light curves covering the 2002 
outburst of 4U 1543-47. 

4.1. WSIM distribution 

Given that each of the representative light curve has 10,000 
points each and given that the WSIM is applied to each 
individual vector, this analysis yields nearly 10,000 values 
of on. The results of this process are illustrated in Figure 3 
which shows the oti distributions for the HS (top panel), 
VHS (middle panel), and LS (bottom panel); the dashed 
lines represent the mean value (a) in the 3 different spectral 
states. 

At first order, all 3 histograms share a similar asym- 
metric shape with a sharp cut-off on the left-hand side and 
a broad right-hand tail. The left-hand side of the on dis- 
tribution is generally related to the correlated variability 
component, whereas the right-hand side is related to the 
random noise component. As a consequence, a highly cor- 
related process is characterized by a narrow on distribution 
peaking at low values. On the other hand, a process dom- 
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Fig. 3. Histograms of the weighted scaling indices (WSIs) 
for HS (top panel), VHS (middle panel), and LS (bottom 
panel). The dashed lines represent the mean values. All 
WSI values were computed for a 3-D embedding space, 
r =0.1, and R = 1.6, using light curves containing 10,000 
data points. 



histogram with a pronounced right-hand tail extending to 
large values of a. 

A first simple way to assess the difference between the 
3 representative histograms is to compare their respective 
means. For the HS, VHS, and LS we get respectively (a) = 
2.038±0.006, 1.769±0.005, 2.013±0.007 (where the quoted 
uncertainties are a/y/n). These values suggest that the VHS 
is significantly different from both the HS (at ~ 33<r level) 
and the LS (28 a level), whereas the difference between HS 
and LS is only marginally significant according to this test 
(2.7a). 

A formal comparison between the 3 distributions based 
on a Kolmogorov-Smirnov test (hereafter K-S test) , which 



the three distributions are statistically different from each 
other. Specifically, the K-S test yields a statistic of 0.22 
and an associated probability Pk-s < 10 -25 that HS and 
VHS on histograms are drawn from the same distribution. 
Similarly, for HS vs. LS and VHS vs. LS we obtain 0.07 
(Pk-s ^ 1CT 24 ) and 0.16 (P K -s < 10~ 25 ), respectively. 
However, it must be kept in mind that the K-S test is 
devised for independent data-points, whereas the different 
on are not completely independent. As a consequence, the 
apparently highly significant difference between the 3 his- 
tograms representing 3 different spectral states should be 
considered with caution and need to be confirmed with fur- 
ther analysis (see below). 

The results from the histogram analysis are encouraging 
and suggest that the WSIM has indeed the potential to dis- 
criminate between the different spectral states of 4U 1543- 
47, in full agreement with results from the PSD analysis. 
However, in order to reach a stronger conclusion, we should 
demonstrate that all oti histograms of HS are indistinguish- 
able from each other, but statistically different from all the 
VHS and LS cti histograms. Although feasible, this proce- 
dure would be very time consuming and would go against 
the primary goal of this work, which is to provide a sim- 
ple alternative way to characterize the temporal properties 
of GBHs. In addition, it must be pointed out that these 
results have been obtained using 10,000-point light curves, 
which are generally not commensurable with typical AGN 
light curves. This approach would therefore hamper a di- 
rect comparison between GBH and AGN, which is one of 
the secondary goals of this work. 

4.2. WSIM mean 

Since our primary goal is to define a simple way to charac- 
terize the global variability properties of GBHs and since in 
this kind of analysis the mean value (a) is the most robust 
indicator of the global variability properties, we will restrict 
our analysis to (a). In this way, the properties of a given 
light curve are defined by a single number, (a), in a similar 
way as the photon index is often used to characterize the 
energy spectral properties of X-ray sources. 

4.2.1. Test with short intervals 

In addition, in order to further generalize this procedure 
extending it to relatively short light curves, which are more 
common than long uninterrupted ones, we will use intervals 
of 100 s (i.e., time series containing 1000 points since the 
bin time is 0.1 s). In this way, the light curves will contain a 
number of points comparable to typical AGN light curves, 
offering the possibility of a direct comparison between GBH 
and AGN variability properties 

Before proceeding further, we must first demonstrate 
that the choice of shorter intervals will not hamper our 
analysis. On one hand, from the PSD analysis of Park et 
al. (2004) we are ensured that nothing relevant occurs to 
the timing properties of 4U 1543-47 for timescales longer 
than ~ 100 s: all the interesting PSD features (frequency 
breaks and QPOs) are located at frequencies well above 
10~ 2 Hz. On the other hand, we must still verify that the 
WSIM results are not significantly affected by the use of 
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Fig. 4. Comparison between the mean scaling index ob- 
tained using long light curves (> 33,000 points; thick pale- 
colored dashed lines) and the values derived using inter- 
vals with 1000 points (smaller dark solid lines), whose av- 
erage is represented by the thick solid lines with horizontal 
lines indicating the dispersion a . The top panel refers to 
a high count rate HS (~16,000 counts/s), whereas the bot- 
tom panel describes results from a low count rate HS (^800 
counts/s). All values of (a) were computed for a 3-D em- 
bedding space, r =0.1, and R = 1.6 
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To this aim, and to demonstrate that the WSIM is 
also independent of the mean count rate (and hence of the 
Poisson noise level), we have performed the following ex- 
periment. We have chosen 2 of the longest HS light curves 
(both have more than 33,000 data points) with mean count 
rate greatly different: the first light curve, obtained close 
to the outburst peak, has an average RXTE PCA count 
rate of ^16,000 counts/s, whereas the second one (corre- 
sponding to a later phase of the decay) has a count rate of 
only ~800 counts/s. We have applied the WSIM to both 
data sets, first using the entire light curve and then using 
intervals of 1000 points only. The results are illustrated in 
Figure [H where the value of (a)33,ooo for the entire light 
curve is represented by the thick dashed line, the individ- 
ual values obtained with 1000 points, (a)ioooj are depicted 
as shorter continuous lines, and their average is indicated 
by the thick continuous line. The horizontal lines indicate 
the standard deviation, <J^ looa , of the sample of (a)iooo- 
Figure |4] reveals that: 

1) In both cases, (a)iooo values narrowly cluster around 
the mean scaling index, (a)33.ooOi obtained from the en- 
tire light curve and their average, ((a)iooo), is fully consis- 
tent with (a) 33,ooo- This is formally demonstrated by the 
fact that |(a) 3 3,ooo~ <<a)iooo)|/(c r ( a ) 1000 /V") = 2.5 (for the 
high count rate case) and 0.1 (for the low count rate case), 
which are both lower than the 3<r level. This indicates that 
WSIM results obtained with 100 s intervals are fully consis- 
tent with those obtained using a longer interval. Therefore, 
with this method the variability properties of 4U 1543-47 
can be thoroughly investigated using 100 s intervals. 

2) Although visually the high and low count rate 
distributions of (a)iooo and their respective mean look 
fairly close and indeed their standard deviations sig- 
nificantly overlap, statistically speaking, their difference 
K(a)iooojiigh) - ((a)iooo,low}| is slightly above the formal 
3cr level. In order to thoroughly address this issue, and es- 
timate quantitatively the uncertainty on the scaling index 
during the HS, we need to account for all the observations 
during HS. This is illustrated in Figure which shows the 
distribution of (a)iooo obtained using all the 100 s inter- 
vals of all the available HS observations. Despite the huge 
difference in count rate (max> 18,000 counts/s, min<370 
counts/s) and a temporal separation longer than 30 days, 
the vast majority of values narrowly clusters around the 
mean, yielding ((a)iooo) = 2.035^^022, where the quoted 
uncertainties are the 90th percentiles (the error on the mean 
is a/s/ri = 0.01/V741 = 0.0004). The "narrowness" of the 
distribution of the (a) values shown in Fig. 5, really sug- 
gests that WSIM results are not affected neither by the fact 
that the light curves have vastly different signal-to-noise ra- 
tios nor by the time span over which the HS light curves 
are spaced. It is therefore suggestive that (a) is most prob- 
ably determined just by one factor, i.e. the properties of 
the variability mechanism during the HS. 



Fig. 5. Distribution of the mean WSI values, obtained us- 
ing segments of 100 s (1000 points) during HS (741 seg- 
ments). All values of (a) were computed for a 3-D embed- 
ding space, r =0.1, and R — 1.6. The thick dashed line 
indicates the mean of the (a) distribution and the dotted 
lines represent the 90th percentiles. 



For completeness, we also carried the previous test for 
the two VHS and LS light curves. Note that the short du- 
ration of the VHS (due to the intrinsic transient behavior 
of this short-lived state) and the LS (due to the interrup- 
tion in the PCA monitoring program) yielded only 2 light 
curves for each state. In both cases, ((a)iooo)) is fully con- 
sistent with (a) whole, as their difference is of the order of 
1 a or less. Similarly, the values of the difference between 
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|((a)iooo,A) - ((Q!)iooo,b)|/V^a" 
the VHS and LS, respectively. 



J, is 0.07 and 1.4 for 



4.2.2. Discrimination of spectral states 

We can now assess whether the WSIM is able to discrim- 
inate between the spectral different states, by considering 
all the available observations, dividing all the individual 
light curves covering the 2002 outburst into 100 s intervals 
and treating each segment as an independent data-set. This 
procedure yields 741 data sets for HS, 31 for VHS, and 40 
for LS. 

Figure O shows the normalized distribution of (a)iooo 
values for the HS (dotted line), LS (solid line), and VHS 
(dashed line), as well as their average of the mean WSIs. 
The HS and LS distributions are very narrow, and they ap- 
pear to be offset, with the LS values of (a) being system- 
atically lower than the respective values of HS. The VHS 
distribution has a rather large width, but is this mainly due 
to the fact that, in addition to the 2 VHS observations, we 
also included here the data from the observation just be- 
fore the source entered fully the VHS state, which caught 
the source during a transition phase (see §4.3). In any case 
though, the VHS values of (a) appear to be systematically 
much lower than the values in the LS and/or HS. 

A simple and robust way to quantify the difference be- 
tween the 3 different states is to compare their respec- 
tive means by computing the quantity Acha-b = K^a) — 
( a B)|/ \J <T \ J r fg, where (a a) is the mean scaling index 
obtained using all the 100 s intervals during the spectral 
state A , and <j\ is the variance divided by the number of 
100 s intervals in that state. This test yields respectively 
AavHS-HS = 8.4, Aa V HS-LS = 7.2, and Aa H s-LS = 10.3, 
indicating that the WSIM is indeed able to statistically dis- 
criminate between the 3 states. 

The significance of the difference between the WSIs in 
the three spectral states can also be examined using a K- 
S test, which can be safely applied to the 3 distributions 
of (oei), since each data-point has been obtained from a 
separate 100 s interval, which in many cases are separated 
by a few days intervals, and hence can be considered as 
being independent "measurements". Specifically, the K-S 
test yields 0.93 (P K s < 8 x 10~ 24 ), 0.78 (P K s < 2 x lO" 10 ), 
and 0.75 (P K s < 6 x 10" 20 ), for the cases of VHS vs. HS, 
VHS vs. LS, and VHS vs. LS. 

In summary, the main results of the analysis based on 
all available light curves divided into 100 s intervals can be 
summarized as follows: 

1) In all 3 spectral states (a) < D (where D = 3 is the 
embedding space dimension). This result implies the pres- 
ence of correlated variability, which is expected given the 
red-noise trend observed in all PSDs. 

2) The 3 spectral states have different mean scaling indices 
satisfying the following relationship: (a)vHS < (o)ls < 
(a)ns- In addition to formally demonstrating that the scal- 
ing index method is able to discriminate between the spec- 
tral states of 4U 1543-47, this results indicates that the 
VHS is the state characterized by the highest degree of 
correlated variability. Also this result is somewhat expected 
from the linear temporal analysis, given the presence of a 
prominent QPO in the VHS PSD. However, the low value of 
(a) vhs might also be related to non-linear correlations, i.e., 




Fig. 6. Normalized distributions and averages of mean 
WSIs for HS, VHS, and LS, obtained using segments of 
100 s (1000 points) from all the available light curves. The 
uncertainties shown represent the error on the mean (i.e., 
a/y/n, where n is the number of data sets); for HS the un- 
certainty is smaller than the thickness of the vertical bar 
representing the mean. All values were computed for a 3-D 
embedding space, r =0.1, and R = 1.6. The maximum of 
the HS histogram has been normalized to 1, whereas the 
VHS and LS have been normalized to 0.15, only for clar- 
ity reasons. No visual comparison would have been possible 
maintaining the real proportions between the distributions, 
given that VHS and LS only contain 30-40 data-points as 
opposed to the 741 contained in the HS 



PSD or autocorrelation function and that manifest them- 
selves as correlation in the Fourier phases. This will be dis- 
cussed in more detail in 55. 



4.3. Temporal evolution of (a) 

After having demonstrated that the WSIM is a reliable 
tool to characterize the variability properties in different 
spectral states, we focus now on the temporal evolution of 
the mean scaling index during the 2002 flare of 4U 1543- 
47. Since the source spends the vast majority of the time 
in the thermal dominated HS, no substantial changes in 
the energy spectral parameters appear to occur during the 
outburst (see Fig. 3 of Park et al. 2004). On one hand, 
this spectral behavior represents a considerable advantage 
to thoroughly assess the uncertainty on (a), but on the 
other hand it partially hampers a detailed study of corre- 
lated spectral and temporal variability. Nevertheless, it is 
interesting to investigate if (and how) the changes of (a) 
are related to the flux changes during the outburst evolu- 
tion. The results of this analysis are illustrated in Figure [TJ 
the top panel describes the temporal evolution of the total 
count rate that is roughly proportional to the evolution of 
disk flux and color temperature; the middle panel shows 
the flux associated with the power law component in units 
of 10~ 10 erg cm~ 2 s _1 , which was derived from Fig. 2 of 
Park et al. (2004); finally, the bottom panel presents the 
(a) temporal evolution. The scaling index analysis is car- 
ried out dividing each individual observation into intervals 
of 100 s, computing (a) for each interval, and then taking 



8 



M. Gliozzi et al.: Nonlinear variability study of 4U 1543-47 




1 B I I . ■ I 

52450 52460 52470 

Time (MJD) 



Fig. 7. Temporal evolution during the 2002 flare of 
4U 1543-47 for the RXTE PCA count rate (top panel), 
flux associated with the PL component in units of 
10 -10 erg cm~ 2 s _1 (middle panel), and mean WSIM (bot- 
tom panel) . The statistical errors are in most cases smaller 
than the symbols. 



The error-bars for the mean scaling indices are given a/y/n 
(where n is the number of 100 s intervals of a given obser- 
vation), which are often smaller than the symbols shown in 
the bottom panel of Fig. [7j 

The temporal trend of (a) can be summarized as fol- 
lows: the mean scaling index remains roughly constant 
around 2.05 for the whole duration of the outburst, with 
the notable exception of a deep dip during the VHS state. 
A closer look at Fig. [Jj reveals 2 minor dips preceding the 
VHS and a steady decrease toward the end of the outburst 
when the source is entering the LS. Interestingly, while 
these changes of (a) appear to be uncorrelated with the 
overall count rate (and hence with the disk flux), they seem 
to occur in correspondence of local maxima in the power 
law flux. At zeroth order, such apparent correlation can be 
interpreted in the following way: the variability associated 
with the power law component is "less random" than the 
one produced by the disk. However, additional data and 
a systematic analysis of correlated spectral and temporal 
properties is necessary before drawing firmer conclusions. 

Perhaps the most remarkable result from Fig. [7] is the 



the total count rate decrease very significantly (in an abrupt 
way the PL flux and in a smoother way the disk flux), (a) 
returns to the same level it was before the VHS state. In 
other words, (a) appears to be a true indicator of state, as 
it is defined based on spectral and timing criteria. 

5. Search for Nonlinearity 

The second goal of this work is to search for nonlinear tem- 
poral correlations in the variability properties of 4U 1543- 
47. Investigating the nature of the temporal variations is 
of primary interest for constraining models of variability in 
GBHs (and AGN). Indeed, unlike the results from the PSD 
and auto-correlation analyses, which can be equally well 
explained by a variety of different physical models, the de- 
tection of nonlinear variability would immediately rule out 
any intrinsically linear model that explain the X-ray vari- 
ability as the superposition of many independent active re- 
gions (e.g., Terrel 1972). Therefore, this analysis has the po- 
tential to provide model-independent constraints that will 
break the current model degeneracy. 

In order to find out whether a time series can be com- 
pletely modeled by superimposed linear processes (plus un- 
correlated noise) or whether signatures of nonlinear correla- 
tions are present, one of the most direct approaches is based 
on the idea of surrogate data sets introduced by Thcilcr ct 
al. (1992; see Kantz, & Schreiber 1997 for a review). In sim- 
ple words, this technique can be summarized as follows: 

1) Assume as null hypothesis that the original time series 
is linear. 

2) Construct linear surrogate data that have the same linear 
characteristics of the original data. In other words, the PSD 
of surrogate data should be indistinguishable from that of 
the original data, since linear processes are by definition 
completely characterized by the PSD or alternatively by 
the autocorrelation function. 

3) Use an appropriate nonlinear statistic to test the null 
hypothesis comparing original and surrogate data. If this 
test yields consistent values for surrogates and real data, 
we conclude that the original time series is linear in nature. 
On the other hand, if the value of the nonlinear statistic 
for the real data is significantly different from the corre- 
sponding values obtained using surrogates, then we infer 
the presence of nonlinearity. 

Before presenting the results from this test, it is instruc- 
tive to describe a few basic details of the technique used to 
create surrogate data and the main characteristics of the 
nonlinear statistic used in this case. 

The most popular algorithms used to generate an 
ensemble of surrogate realizations are the Amplitude 
Adjusted Fourier Transform (AAFT) and the Iterative 
Amplitude Adjusted Fourier Transform (IAAFT) algo- 
rithms (Theiler et al. 1992; Schreiber & Schmitz 1996). 
Although the AAFT and IAAFT algorithms conserve the 
amplitude distribution in real space and reproduce the PSD 
of the original data set quite accurately, it has been shown 
recently that both algorithms may induce unwanted cor- 
relations in the Fourier phases (Rath & Monetti 2008). To 
guarantee that the surrogates in this study are free from any 
higher order correlations, we generate them in the follow- 
ing way: First, the time series is mapped onto a Gaussian 
distribution in a rank-ordered way, which means that the 
amplitude distribution of the original times series in real 
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Fig. 8. Nonlinear prediction error (NLPE) for the HS (top 
panel), VHS (middle panel), and LS (bottom panel). The 
black lines show the NLPE for the original time series. The 
smaller thin colored lines denote ip for the respective surro- 
gates. The thicker and longer dashed line indicate the mean 
value (ip) as derived from the 50 surrogate realizations. The 
dash-dotted lines mark the 3<7^, interval, i.e. (ip) — 3a,p and 
(ip) + 3(7.0 , where the standard deviation is also derived 
from the surrogates. 



that the rank-ordering is preserved, i.e. the lowest value of 
the original distribution is replaced with the lowest value of 
the Gaussian distribution etc. By applying this remapping 
we automatically focus on the temporal correlations of the 
data, while excluding any contributions to nonlinear cor- 
relations stemming from the non-Gaussianity of the origi- 
nal intensity distribution. Second, we Fourier transform the 
remapped time series, replace the original phases by a new 
set of uniformly distributed and uncorrelated phases and 
perform an inverse Fourier transformation. Note that the 
surrogate time series generated in such a way preserve ex- 
actly the power spectrum, while explicitly controlling the 
randomness of the phases. For each time series under study 
we generate 50 corresponding surrogates. 

Although in principle any nonlinear statistics may be 



performed by Schreiber & Schmitz (1997) indicates that 
the nonlinear prediction error (NLPE) is one of the most 
effective indicators of nonlinearity and has good discrimina- 
tion power to detect any deviation from a Gaussian linear 
stochastic process. Therefore, to test the presence of non- 
linearity in 4U 1543-47 light curves we make use of the 
NLPE. 

The equation used to compute the NLPE (hereafter ip) 
as well as some technical details are provided in Appendix 
B. Here we simply describe in a qualitative way the 
main characteristics of this nonlinear indicator. Since this 
method relies on the time delay embedding technique de- 
scribed in §3, it can exploit the direct correspondence be- 
tween the scalars of the starting time series, x\, ■■■,x n , 
and the corresponding set of vectors in the pseudo phase 
space: Xi, ,x 2 , ...,x n . Specifically, to predict the "future" 
measurement x n +io (i.e., the value of the time series 10 
steps ahead of x n ), one must find the closest vector to x n , 
which we will call Xi and corresponds to the scalar Xi, and 
then use the scalar Xi + io as a predictor for x n+ iQ. As ex- 
plained in Appendix B, the rigorous process is slightly more 
complicated and involves several vectors in the neighbor- 
hood of x n , which in turn will yield several predictors. The 
nonlinear prediction error ip is then provided by the aver- 
age of the differences between the actual value x ra +io and 
the different predictor values. 

We have carried out this test for all the light curves 
relative to the 2002 outburst of 4U 1543-47, although for 
clarity reasons in Figure 8 we only show the results for the 
3 representative light curves of the HS, VHS, and LS intro- 
duced in §3. From this figure, it is clear that the absolute 
values of ip for the HS and IS are considerably larger than 
those measured in the VHS. This is an expected behavior, 
because in the VHS the time series is much more correlated 
with linear components showing up as QPO. As a conse- 
quence, its predictability increases leading to lower values of 
ip. More importantly, Figure 8 reveals that the VHS shows 
highly significant signatures for nonlinear correlations, un- 
like the HS and IS for which the values of ip obtained using 
real data are fully consistent with the values derived using 
linear surrogates. 

Analyzing all time series belonging to the different spec- 
tral states we infer that all the HS and LS light curves are 
linear, whereas the 2 VHS light curves show respectively 
highly significant and marginally significant signs of non- 
linearity. Specifically, the time series on MJD=52461, which 
corresponds to the absolute minimum of the scaling index, 
shows the strongest and the only statistically significant evi- 
dence (~ 7.5c) for nonlinearity, as illustrated by the middle 
panel of Figure 8. A day before, on MJD=52460, the evi- 
dence of nonlinearity is marginal (at a ~ 2.4c significance 
level). Finally, the light curve on MJD=52459 that caught 
4U 1543-47 during the HS-to-VHS transition (see Fig. 7) 
shows no indications of nonlinearity at all. 

In conclusion, the surrogate test reveals that: 1) All HS 
and LS light curves are linear, 2) non-linearity indications 
appear during the VHS light curves, and the highest and 
most significant signal for non-linearity occurs for this light 
curves with the strongest QPO and the lowest value of (a). 
On one hand, the latter result suggests that the physical 
mechanism leading to strong QPOs is intrinsically nonlin- 
ear, implicitly disfavoring QPO linear models. On the other 
hand, it indicates that the low value of (a) measured in the 
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ear correlations, which cannot be detected by linear timing 
techniques. 

6. Summary and Conclusions 

We have carried out a nonlinear analysis of the variability 
properties of the X-ray nova GBH 4U 1543-47. The main 
results can be summarized as follows: 

— We have used the WSIM to assign a single number, the 
mean scaling index (a) , to each individual light curve of 
the source. The large number of data when the source 
was in its HS, showed that the resulting (a) values re- 
main roughly constant, irrespective of large changes in 
flux associated to the disk and PL components. 

— Similarly, the mean scaling index values remained 
roughly constant during the VHS and LS, showing the 
following relationship: (a)vHS < («)ls < (o:)hs- 
These results, which need to be confirmed for other 
GBHs, suggest that the mean scaling index (a) may 
be used to parametrize the timing properties of an ac- 
creting source, and that it may be a true indicator of 
"state" in these systems. 

— When plotted versus time, the (a) trend shows no direct 
correlation with the total flux temporal behavior, which 
is dominated by the accretion disk emission. On the 
other hand, the temporal evolution of (a) appears to 
be somewhat related to that of the power-law spectral 
component: (a) reaches its absolute minimum roughly 
at the same time as the PL flux reaches its absolute 
maximum, which occurs during the VHS state. 

— The search for nonlinearity using surrogate data and 
NLPE reveals that the variability is linear in all light 
curves with the notable exception of one observation in 
the VHS, which corresponds to the absolute minimum 
of (a) and is also characterized by the presence of a 
strong QPO. 

An important implication of the detection of nonlinear- 
ity in the VHS is that all intrinsically linear models pro- 
posed to produced low frequency QPOs (LFQPOs) may be 
ruled out for 4U 1543-47. The formation of QPOs (at both 
high and low frequencies) in GBHs is currently a matter of 
debate, as several viable competing models have been pro- 
posed (see, e.g., McClintock & Remillard 2006 for a review). 
It is therefore important to derive model-independent con- 
straints that may break or at least reduce the current model 
degeneracy. If the detection of nonlinearity associated with 
strong LFQPOs is confirmed in other GBHs, then intrin- 
sically linear models such as the global disk oscillations 
proposed by Titarchuk & Osherovich (2000) can be safely 
ruled out, whereas alternative models such as the shock 
oscillation model (Chakrabarti & Manickam 2000) or the 
accretion-ejection instability model (Tagger & Pellat 1999), 
which may account for nonlinearities, are still viable solu- 
tions. 

In summary, the findings derived from this work sug- 
gest that this kind of nonlinear analysis can be useful in 
the field of GBHs and complement the temporal analysis 
carried out with linear techniques. Specifically, the simplic- 
ity of the WSIM, which characterizes the global variability 
properties of a light curve via a single number, (a) , suggests 
that this technique might be successfully applied to study 



the robustness of the WSIM, that performs well also with 
noisy data and relatively short light curves, naturally sug- 
gests an application to AGN variability with the possibility 
of direct comparison with GBHs. 

Since we have limited our analysis to 4U 1543-47 only, 
before deriving any general conclusions, we should carry 
out a systematic nonlinear analysis on several GBHs during 
their spectral transition. Specifically, it will be important 
to assess whether specific values of (a) are associated with 
specific spectral states (e.g., (a)vHS — L8, (a)us — 2.0, 
(a)ns — 2.05), or if only the temporal trend of (a) during 
the outburst (i.e., the presence of a pronounced minimum 
during the VHS) is similar to the one displayed by 4U 1543- 
47. To this aim we will carry out a similar analysis on a few 
prominent GBHs with outbursts that are completely cov- 
ered by the RXTE PCA and where the different spectral 
state are well sampled. This test will unequivocally reveal 
whether (a) can be used as reliable indicator of spectral 
states. In addition, the planned study will provide useful 
information on correlated variability of (a) and several rel- 
evant spectral parameters, and on the presence of nonlin- 
earity in different spectral states. 
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Appendix A: Dependence of SIM on r, D, and R 

As explained in §3, the WSIM depends on 3 parameters: 
time delay r and embedding dimension D, which are related 
to the phase space reconstruction process, and the radius R 
at which the logarithmic derivative (i.e., the scaling index) 
is computed. 

In this section, we visually demonstrate that the impact 
of these parameters on our main findings is basically irrel- 
evant for a broad range of reasonable choices. Figure. IA.1I 
shows the temporal evolution of (a) for D =2 (top panel), 
3 (middle panel), and 4 (bottom panel). Solid lines refer 
to r = 0.1 s, dotted lines to r = 1 s, and dashed lines to 
t = 100 s; in this case the radius is fixed at 1.6. It is no- 
ticeable that all the values of (a) consistently increase as 
the embedding dimension increases. This reflects the fact 
that a significant part of the variability is random and this 
random component translates into larger values of a as the 
dimension of the embedding space increases. 

From Fig. IA.1I it is evident that all different combina- 
tions of r and D reproduce the same temporal evolution 
trend of (a) with a large central dip preceded by to low- 
amplitude dips and followed by a final steady decrease. It 
is also clear that using long time delays, such as r = 100 
s, the dips appear less prominent. This simply reflects the 
fact that long time delays necessarily loose the information 
relative to short-term temporal correlations. 

Figure. [A. 2l illustrates the impact of the radius R on the 
temporal evolution of (a) (this time D = 3 and r = 0.1 s). 
Once again, all values of R are able to recover the same 
temporal evolution trend of (a) described above. 

Appendix B: Nonlinear Predictor Error 

To calculate NLPE, the time series is embedded in a D- 
dimensional space using the method of delay coordinates 
as described in section 3. We use here also the embedding 
dimension D = 3. The delay time r was determined using 
the criterion of zero crossing of the autocorrelation function 
considering the VHS, where the ACF is sufficiently different 
from a random process. The NLPE is defined as 

M-l-T \ 

£ [*„ +T -F(*„)H (B 

n=(D-l)r / 

where F is a locally constant predictor (i.e., a quantity that 
remains constant for a local surrounding of a point under 
study in the D-dimensional embedding space), M is the 
length of the time series, and T is the lead time (i.e., the 
number of time steps ahead of the considered one for which 
we want to make a prediction). The predictor F is calcu- 
lated by averaging over future values of the N (N = D + 1) 
nearest neighbors in the delay coordinate representation. 
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Fig. A.l. Temporal evolution of the mean WSIM during 
the 2002 flare of 4U 1543-47 for embedding dimensions 2 
(top panel), 3 (middle panel), and 4 (bottom panel). Solid, 
dotted, and dashed lines refer to r = 0.1, 1, and 10 s, 
respectively. 
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Fig. A. 2. Temporal evolution of the mean WSIM during 
the 2002 flare of 4U 1543-47 for embedding dimensions 3, 
with r ranging from 0.6 (top dotted line) to 2.4 (bottom 
dashed line). 
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time T and did not find significant variations. We set this 
value to T — 10 time steps. The results of this analysis for 
the 3 representative light curves of the HS, VHS, and LS 
are shown in Figure 8. 



